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Epidemics seldom occur as isolated phenomena. Typically, two or more viral agents spread within the same 
host population and may interact dynamically with each other. We present a general model where two viral 
agents interact via an immunity mechanism as they propagate simultaneously on two networks connecting the 
same set of nodes. Exploiting a correspondence between the propagation dynamics and a dynamical process 
performing progressive network generation, we develop an analytic approach that accurately captures the dy- 
namical interaction between epidemics on overlay networks. The formalism allows for overlay networks with 
arbitrary joint degree distribution and overlap. To illustrate the versatility of our approach, we consider a hy- 
pothetical delayed intervention scenario in which an immunizing agent is disseminated in a host population to 
hinder the propagation of an undesirable agent (e.g. the spread of preventive information in the context of an 
emerging infectious disease). 

PACS numbers: 87.23.Ge, 89.75.Hc, 87.10.Ed 



I. INTRODUCTION 

Epidemic dynamics has been largely studied with the help 
of mathematical models in which a single viral agent propa- 
gates in a given host population. Although the paradigm of 
isolated epidemics may be well suited in some cases, there 
are numerous other situations in which more than one process 
occurs and interacts in the same population. Different bio- 
logical pathogens may interact through ecological |Q1 0] and 
immunological Q, [|] mechanisms, or multiple strains of the 
same disease may compete for hosts according to some cross- 
immunity profile lUSl- The spread of fear or awareness in 
the context of an emerging disease |p7l— TT"oh can also be consid- 
ered as a case of interacting viral agents, i.e., information and 
disease. In computer networks, the dissemination of coun- 
termeasures using a contagious vaccination scheme has been 
suggested to counter harmful computer viruses llilT - 11311 . 

When propagating in some host population, two viral 
agents may follow different - or share similar - routes of 
transmission. Taking into account how individuals are in con- 
tact with each other then becomes of great importance when 
modeling their interaction. By explicitly considering those 
heterogeneous contact patterns between individuals, network- 
based models are an ideal framework for the study of interact- 
ing epidemics in structured populations ITil [T5ll . 

The interaction between two viral agents has been studied 
from the perspective of complex networks in a limited number 
of contributions J^-HH [l6l - l20ll . An important step towards a 
general theory of interacting processes on complex networks 
was recently accomplished by Funk and Jansen ll20ll . Gener- 
alizing the previous work of Newman |[l6ll . they analyzed the 
interaction between two viral agents propagating successively 
on overlay networks, i.e., two networks connecting the same 
set of nodes. Albeit very elegant, their analytical approach, 
based on bond percolation, is static and does not apply to the 
case of dynamically interacting viral agents. 



* http://dynamica.phy.ulaval.ca 



The purpose of this contribution is to develop an analyt- 
ical approach able to capture the dynamical interaction be- 
tween viral agents spreading simultaneously on overlay net- 
works. To this end, we make use of a correspondence between 
propagation on networks and a dynamical process perform- 
ing progressive network generation ll2lil (see also Appendix C 
of I22IO . The formalism obtained is quite general, and allows 
for overlay networks with arbitrary joint degree distribution 
and overlap. The language of epidemiology is used through- 
out this work for its clarity, yet our approach may be applied 
to processes of other nature that spread on networks. 

This paper is organized as follows. In Sec. [II] we intro- 
duce a model in which two viral agents propagate and interact 
dynamically on overlay networks. The analytical approach is 
then developed in two steps in Sec. [Ill] Some properties of 
the model are investigated in Sec. [IV] where we also validate 
the accuracy of the analytical predictions by comparison with 
Monte Carlo simulations of the dynamics. Our conclusions 
are summarized in Sec. [V] An online supplementary docu- 
ment J23ll containing the full set of equations of our analytical 
approach completes this contribution. 



II. INTERACTING EPIDEMICS ON OVERLAY 
NETWORKS 

In this section, we introduce a general epidemic model 
where two dynamically interacting viral agents spread on 
overlay networks. 



A. Overlay networks 

A system of two overlay networks is defined by two net- 
works, Y\ — ( V, Ei) and F2 = (V, E%), which connect the same 
set of nodes V through their own set of links, E\ and Eo [see 
Fig. |l(a)| I20I1 . Nodes represent individuals of a given host 
population, while links correpond to potential transmission 
routes between pairs of individuals. As usually done in the 
network litterature, we denote the size of the host population 
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FIG. 1. (a) Schematical illustration of a system of two overlay net- 
works, T[ = (V, Ei) and T 2 = (V, E 2 ). (b) This system can be de- 
composed into three non-overlapping networks: y { = (V, Ei \ E 2 ), 
y 2 = (Y,E 2 \ E l ), and y„ = (V, E l n E 2 ). 

by N = | VI. Two nodes are said to be neighbors on T g , where 
g e {1,2}, if they share a link on this network, and the num- 
ber of neighbors of a node on F g define its degree k g on this 
network. Neither can a node be linked to itself (no self-loops) 
nor share more than one link with another node on the same 
network (no repeated links). 

The system is characterized by its joint degree distribution 
P{kx , ki) = P(k), which corresponds to the probability that a 
node selected at random in the host population has a degree k\ 
on Fi and on Ft. The marginal degree distribution P g (k g ) 
of each network can be obtained by summing over P(k), i.e., 

Piiki) = YjP(h,h) , P2(k 2 ) = Y J P(kuk 2 ). (1) 

k 2 ft 

Even if E\ + E2 in the general case, a given number of links 
may be common to both networks. Because P(k) contains no 
information about the potential overlap in the system, we will 
resort to a useful decomposition into three non-overlapping 
networks lEoll . Let y\ — (V, E\ \ E2) be the network char- 
acterized by all the links unique to Fi, jj_ = (V, E2 \ £1) be 
the network characterized by all the links unique to F2, and 
jb - (V, Ei n En) be the network containing all links common 
to Fi and F2 [see Fig. 1 1 (b)| . Because it contains information 
about the distribution of overlapping links, the joint degree 
distribution resulting from the three-networks decomposition, 
denoted by p{K\ , K2, Kb) = p(k), offers a more accurate descrip- 
tion of the system. The distribution P(k) can be obtained from 
P(k) by summing over Kb, 

P(k) = ^ p(ki -K b ,k 2 -K b , K b ) . (2) 

Kb 

B. Interacting epidemics 

In the single viral agent susceptible-infectious-recovered 
(SIR) dynamics on a network, nodes are divided into three 
states: susceptible (S), infectious (/), or recovered (R). Infec- 
tious contacts occur between infectious nodes and their neigh- 
bors at the rate (3. If a susceptible node is involved in an infec- 
tious contact, transmission ensues and it becomes infectious. 
Infectious nodes recover at the rate a, and become immune to 
further infection. 

Here, we extend the SIR dynamics to two interacting viral 
agents. We assume that "agent 1" propagates via the links 



of Fi, while "agent 2" spreads on F2. If both viral agents 
are transmitted through the same type of interactions between 
individuals, then Fi = F2; otherwise, Fi + 

At any time, the state of a node is given by the combination 
of its particular state regarding each viral agent. For agent g, 
its g-state can be either g-susceptible (S g ), g-infectious (I g ), 
or g-recovered (R g ). The rate of infectious contacts on F g is 
fl g , while g-infectious nodes recover at the rate a g and become 
immune to further infection by agent g. In order to study the 
dynamical interaction between epidemics, we consider a case 
of leaky partial immunity flUl . When an infectious contact 
occurs between a g-infectious node and a g-susceptible node 
whose state regarding the other agent g is Y g e {S g ,I g ,R g }, 
transmission of agent g successfully follows with probability 
crj and the g-susceptible node becomes g-infectious. Other- 
wise, it remains g-susceptible with complementary probabil- 
ity <xj = 1 - crj. We assume that an infectious contact on T g 
may occur only once between two given nodes. 

The motivation for this interaction rule is that, while re- 
maining general, it renders the model analytically tractable us- 
ing a reasonable level of complexity. Other interaction mech- 
anisms could have been considered, such as perfect partial im- 
munity BI81I or leaky partial immunity allowing for more than 
one infectious contact between two nodes. More on this topic 
is covered in Sec. MI CI 



C. Monte Carlo simulations 

Monte Carlo simulations of epidemic propagation on over- 
lay networks are performed in two steps, network generation 
and viral agent propagation. Our overlay networks are gen- 
erated using two different algorithms, both adapted from the 
well-known configuration model ll24ll25ll . 

The first algorithm, based on the joint degree distribution 
P(k), generates two overlay networks with random overlap. 
(i) A random degree sequence {£,} of length N subjected to 
P(k) is generated. Since a link consists of two stubs, we en- 
sure that 2i^i,g i s even for all g e {1,2}, otherwise an ele- 
ment of the degree sequence is selected at random and gener- 
ated again; (ii) For each ki, a node with stubs on Ti and 
£, 2 stubs on T2 is created; (iii) Independently for each net- 
work, pairs of unconnected stubs are randomly chosen and 
connected together until all unconnected stubs are exhausted; 
(iv) The presence of self-loops and repeated links is tested on 
each network. All faulty links on a network are removed by 
randomly choosing a pair of connected stubs on the same net- 
work and rewiring them to the former stubs. 

To generate networks with arbitrary overlap, we use a sec- 
ond algorithm based on the joint degree distribution p{K). The 
procedure is essentially the same as described above, except 
that three networks {y\, yo, and jb) are generated from p{K) 
with the additional constraint that one link cannot exist on 
more than one network. The networks Fi and F2 are then con- 
structed from ji, j2, and 77, as illustrated in Fig.Q] 

The viral agent propagation phase is carried out using dis- 
crete time steps of length At. At each time step and for viral 
agent g e {1,2}, every link on F g between a g-infectious and a 
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g-susceptible node is tested for infectious contacts with proba- 
bility /? ? Af, on the condition that it did not happen previously. 
If the test returns positive, the g-susceptible node becomes 
^-infectious with transmission probability cr Y , where Y g cor- 
responds to its state regarding the other agent g. Recovery 
events are tested with probability a g At. All simulations are 
initialized by infecting at random a fraction e\ of the nodes in 
the system with agent 1 at time t — 0. To allow for a delay 
between epidemics, the dynamics of agent 2 are initialized at 
time t — t > with the random infection of a fraction ei of 
the nodes in the system. 

The simulations presented in this paper are carried out on 
networks of size N = 25 000 (unless explicitly noted) with 
At = 0.001. Both recovery rates are set to unity, a\ — a2 = 1. 
Moreover, we will use e\ — 62 = 0.001 as initial conditions. 
All Monte Carlo results shown in the figures are computed 
over a total of 1000 simulations unless explicitly noted. 



III. MEAN-FIELD APPROACH 

We now develop a network-based compartmental formal- 
ism that describes the dynamics of the model introduced in 
the previous section. Our approach is based on the concept 
of on the fly network generation recently introduced in lull . 
The case of interacting epidemics on overlay networks with 
random overlap is treated in the first place, and the approach 
is later generalized to allow for arbitrary overlap. Finally, we 
provide some insights on how to handle similar types of inter- 
action dynamics. 



A. Networks with random overlap 

In order to develop the analytical approach that follows, 
some preliminary considerations are needed. In Sec. Ill CI we 
explained the two-step procedure used in Monte Carlo simu- 
lations of epidemics on configuration model networks. How- 
ever, as we recently pointed out in a recent contribution 1I2TI1 . 
an alternative procedure, consisting of one single step where 
the networks are generated on the fly (i.e., during propagation, 
when required), can be considered. 

Recall that in the procedure of Sec. Ill CI we test for infec- 
tious contact every link on T g between a ^-susceptible and a 
g-infectious node, on the condition that an infectious contact 
never happened between them. Suppose instead we used to 
test for infectious contact every stub on T g that emanates from 
a ^-infectious node and that has never been the host of an in- 
fectious contact, without any further distinction. Then only 
when this test returned positive would we have inquired about 
the state of the corresponding neighbor. If g-susceptible, then 
we would have tested for transmission. Unlike in Sec. Ill CI 
this new procedure does not require any explicit prior knowl- 
edge of the structure of the networks in the system: neighbors 
can be assigned on the fly by matching stubs pairwise at the 
moment infectious contacts occur. A schematization of on the 
fly network generation for a single agent SIR epidemic is il- 
lustrated in Fig. 12 
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FIG. 2. Schematical illustration of on the fly network generation 
for a single viral agent susceptible-infectious-removed epidemic (Q: 
susceptible nodes, ®: infectious nodes, •: recovered nodes), (a) 
Initially, one node of degree 2 is infectious, (b) When an infectious 
contact occurs over one stub belonging to the infectious node, it is 
matched with another stub chosen at random between all the un- 
matched stubs (dashed lines), thus forming a link (solid lines). If 
previously susceptible, the assigned neighbor becomes infectious, 
(c) Stubs are progressively matched pairwise as infectious contacts 
occur in the population. In the matching process, self-loops and re- 
peated links are allowed, but their probability decreases as yV~'. (d) 
The epidemic stops when there are no more infectious nodes. Parts 
of the network not reached by the viral agent are never built. 



In the case of overlay networks with random overlap, i.e., 
when stubs are matched independently at random on T\ and 
F2, an almost exact correspondence may be established be- 
tween the algorithm exposed in Sec. Ill CI and an equivalent 
stochastic Markov process performing on the fly network gen- 
eration l2lll . The only difference arises from self-loops and 
repeated links, which are allowed in the stochastic process. 
However, since their probability decreases as N , we expect 
the results to agree in the limit of large networks. 

Instead of tracking the full stochastic process, we rather fo- 
cus on mean values in the asymptotic limit N —> 00 and hence 
obtain a fully deterministic approach. Let {XiY^ijif) be the 
mean fraction of nodes in the system that are of 1 -state X\, 
of 2-state Yi, have i unmatched stubs on Fj, and j unmatched 
stubs on T2 at time t l26ll . The ordinary differential equa- 
tion (ODE) governing the time evolution of the [XiY^ly com- 
partment consists of two parts, accounting for the dynamics 
of each viral agent. As both parts are very similar for cor- 
responding states (e.g., agent 1 dynamics for 1 -susceptible 
nodes and agent 2 dynamics for 2-susceptible nodes), we ex- 
clusively focus on agent 1 dynamics. The same considerations 
apply to agent 2 dynamics as well. 

Let 0i be the probability that an unmatched stub on Ti be- 
longs to a 1 -infectious node, 



0, = 



(3) 



Nodes in the [Si^Ly compartment will lose unmatched 
stubs on F! at the rate /3\®\i as they are involved in infec- 
tious contacts with 1 -infectious individuals. When this hap- 
pens, they are either transferred to the [/iF2](/-i)/ compart- 
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ment with probability cr\ (successful transmission), or to the 
[S i Y2\t-i)j compartment with probability <xf (failed transmis- 
sion). This yields the following contribution to the ODE gov- 
erning [SiY 2 ]ij: 



0i0i[(»+l)af[SiY2]<WM- i[SiY 2 ] u ] 



(4) 



Nodes are transferred from [/iF 2 ];/ to [i?iF 2 ]y at the rate ffi 
due to recovery events. Moreover, a node from [I\Y2\ij can be 
transferred to [/il2](i-i)j due to a loss of an unmatched stub, 
which may occur in two different ways: if it is the source 
of an infectious contact (at rate or if it is victim of an 
infectious contact originating from another 1 -infectious node 
(at rate j3\®\i). Adding the incoming flow of newly infectious 
1 -susceptible nodes, one obtains the following contribution of 
agent 1 dynamics to the ODE governing [/liy*,-: 



-ai[/iF 2 ]y+A®ia+l)o-[[SiF 2 ] (i+1)j 
+P l (l+& l )[(i+l)[I l Y 2 ] (i+l)j - i[/lF 2 ]y] . 



(5) 



Nodes are removed from [/^i^Ly an d transferred to 
[R\Y2](i-\y at the rate [5\<d\i as they are involved in infectious 
contacts with 1 -infectious individuals. Including the incoming 
flow of newly recovered 1 -infectious nodes, the contribution 
of agent 1 dynamics to the ODE governing |7?iT 2 ],y reads 

ait/i^y+jSiQi^'+^^i^ci+D-^i^y] . (6) 

Finally, in order for the dynamics to be completely spec- 
ified, an initial condition is required for each compartment. 
As mentioned earlier, random infection of a fraction e\ of the 
population with agent 1 occurs at t — 0, which gives: 



eiP(i,J) 




if X = S and Y = S 
if X = /and Y = S 
otherwise . 



(7) 



Agent 2 is then introduced at random in the population at time 
t — t and infects a fraction e 2 of the nodes. This is imple- 
mented by making the following substitution at t — r: 



[XiY 2 Mt) 



(l-«*)[XiS 2 ]y(T) ify=S 
^[XiSdyW ifF=/ (8) 
otherwise . 



The full system of ODEs describing the dynamics of the 
model on overlay networks with random overlap can be found 
in the supplementary document of this paper 023I1 . 



B. Networks with arbitrary overlap 

Generalization to the case of arbitrary overlap requires ad- 
ditional considerations. They originate from the fact that the 
two networks are not built independently anymore by match- 
ing stubs at random on each network. Since the structure of 
one network is now influenced by the structure of the other, 



an exact correspondence cannot be established with the dy- 
namical process performing on the fly neighbor assignment 
introduced in Sec. MI Al . We therefore have to rely on some 
approximations. 

Consider the networks j\, y 2 , and % resulting from the 
three-networks decomposition presented in Sec. Ill Al Let 
[X\Yz]ijk(t) be the mean fraction of nodes in the system that 
are of 1 -state X\, of 2-state F 2 , have ; unmatched stubs on y\, 
j unmatched stubs on y 2 , and k unmatched stubs on yi, at time 
t. Once again, let us concentrate on the part of each ODE that 
corresponds to agent 1 dynamics. 

In Sec. 1111 Al the case of random overlap was considered 
in the asymptotic limit N — > oo. Since one has Y\ = y\ and 
F 2 = j2 under this condition, substituting [XiF 2 ], ; by [-XiF 2 ]yj 
in Eqs. ©-(|6]l yields the correct description of agent 1 dy- 
namics on j\ . Hence, we only need to derive the additional 
contributions corresponding to the dynamics of agent 1 on yt,. 

The approximation that we will use to take into account the 
overlap between Y\ and T 2 can be illustrated by the follow- 
ing example. Consider two nodes, node n (of state \X\Y2\ji) 
and node n' (of state [X' { YL]tfkf)> that are neighbors on H and 
T 2 , and that are just being involved together in an infectious 
contact with agent 1. From the point of view of a dynami- 
cal process performing on the fly neighbor assignment, their 
respective number k and k' of unmatched stubs on yi, should 
be decreased by one. However, n and n' could be later in- 
volved together in an agent 2 infectious contact. Since we do 
not track the information about the states of a node's neigh- 
bors in the formalism, one could account for this by increas- 
ing their number j and / of unmatched stubs on y 2 by one 
at the same moment k and k' are decreased. In other words, 
the prior information that n and n' were neighbors would be 
forgotten, but the fact that they may later be involved together 
in an agent 2 infectious contact is approximatively accounted 
for by granting them a new unmatched stub on y 2 . If an in- 
fectious contact occurs later over one of those stubs, it will be 
matched with another stub chosen at random between all the 
unmatched stubs on y 2 . 

Yet by doing so, we are throwing some useful information 
away. When two nodes are involved in an infectious con- 
tact, they share information about their respective states. If 
node n was [/lSJy* an d node n' was [S \R 2 ]i>fk', the knowl- 
edge that node n' is 2-recovered tells us that a transmis- 
sion of agent 2 will never occur later between them. Giv- 
ing them an unmatched stub on y 2 would thus be ill-advised; 
the ijk — > ij(k-l) and i'j'k' — » i'j'(k'-l) transitions would 
be more appropriate. In summary, if node n is involved in 
an agent 1 infectious contact via y\, with node «', they will 
be given an additional unmatched stub on y 2 if and only if 
their states are such that an agent 2 transmission may occur 
between them later in time. 

Let Q Y bx be the probability that an unmatched stub on y\, 
belongs to a 1 -infectious node of 2-state F 2 : 



£1 



Yjij,kk[hY2]ijk 



(9) 



The total probability 0/,j that an unmatched stub on yt, be- 
longs to a 1 -infectious node is then given by ©/, 1 = £)y 6? v 



5 



We further define the probability ©J l that an unmatched stub 
on yi, belongs to a 1 -infectious node whose 2-state is such 
that a transmission of agent 2 may occur later with a node of 
2-state Y2: 







I, A 



if Y-- 
if Y-- 



(10) 



— v 

The probability of the opposite event is & h x = 0/,,i — 0f v 

Similarly to Eq. ||9), we define the probability <f>Y j that any 
unmatched stub on y\, belongs to a node of particular 2-state 
Y 2 : 



Zx Yjij,kk[XiY 2 ]ijk 



(11) 



Naturally, 'Zjy'PIi = 1- The probability OT^ that an un- 
matched stub on yb belongs to any node whose 2-state is such 
that a transmission of agent 2 may occur later with a node of 
2-state Y2 is given by: 



/>.! 



b' bl if Y=S 
if Y=I 
\fY=R 



(12) 



We denote the complementary probability <& hl = \ — OL. 

First, a node in the [5 li^lrjA compartment loses unmatched 
stubs on yi, at the rate /3\®b,\k as it suffers infectious con- 
tacts with 1 -infectious individuals. This rate is composed 
of two distinct parts: infectious contacts from 1 -infectious 
nodes with whom a transmission of agent 2 may occur later 
{f3\® Y h] k) and from 1-infectious nodes with whom a transmis- 

— Y 

sion of agent 2 may never occur (J3\Q b l k). As mentioned ear- 
lier, the first part yields a ijk — * i(J+Y)(k—Y) index transition, 
while the second gives the transition ijk — » ij(k- 1). After 
an infectious contact, the [S1I2];;* node becomes 1-infectious 
with probability crY, or remains 1 -susceptible with probability 
erf. Combining those four different issues, the contribution of 
agent 1 dynamics on yt, to the ODE governing [S 1 F2L7/S: is 



-p x ® bA k[S x Yi\uk 



y-l)(i+l) + ®b,l$ l Y 2\ij(k+l)J 



(13) 



Second, a node in the [/1 ^2]/^ compartment loses un- 
matched stubs on yb at the rate /3i(l +®b,i)k as it is the source 
or target of an infectious contact. Once again, this rate is com- 
posed of two parts that yield different index transitions: infec- 
tious contacts with a node with whom a transmission of agent 

2 may [fiii®^ + ® Y bl )k] or may not + ®l A )k] occur 

later. Taking into account the flows of incoming 1 -susceptible 
nodes, the contribution of agent 1 dynamics on yb to ODE for 



(14) 



UiYilijk reads 

p x {k+\)oi[®llS , Y 2 ] iU - m+ i } + 0^SiF 2 ]jAH-i)) 

+Py(k+\)i^llhY2]i { j- m+ l)+^lihY2} im ^ 

+j8i(*+l) (^f/ift]iCMX)h-i) + ©m[A ^1,7(^1)) 
-Ad+Oi.i^/iTz],^. 



Finally, a node in the [R\ Y%\ijx compartment loses unmatched 
stubs on yb at the rate /3\®b,\k as it suffers infectious contacts 
with 1-infectious individuals. Taking into account both possi- 
ble issues, i.e., if an unmatched stub is added afterwards on 72 
or not, one obtains the contribution of agent 1 dynamics on yb 
to the ODE governing [R\Yz\ijk- 



—Y 

U-w+i) + &b y i[RiY 2 ]ij(k 



-PiQbMWlijk 



(15) 



The interested reader is once again referred to the supplemen- 
tary document for the complete system of ODEs 112311 . 

Initial conditions for agent 1 dynamics are given at t = by 



Ry 2 ]i#) 



(l-£i)p(i,j,k) ifX = SandT = S 
eip(i,j,k) ifX = /andT = 5 

otherwise , 



(16) 



while agent 2 dynamics is initialized at t — r with 

[XiY 2 ] U k(T) 



(i-ft)[XiS 2 ] j/t (T) if y=s 

e2[X 1 S 2 ] i jk(T) if Y=I 
otherwise 



(17) 



The complexity of the ODE system derived in the case of 
random overlap increases as <9(A: l | nax xA:™ ax ), while it increases 
as <9(/c™ x XA-™ ax x/c™ ax ) in the case of arbitrary overlap. Here 
A: max and *r max represent the largest degrees beyond which the 
systems of ODEs are truncated. Although high complexity 
may seem to be a major drawback of our approach, a signifi- 
cant speed up on Monte Carlo simulations can be obtained. 



C. Other types of dynamics 

As we have pointed out at the end of Sec. Ill Bt other types 
of interaction dynamics could have been considered. We give 
two examples to show how such alternative interaction rules 
could have been analytically handled using similar modeling 
schemes. 

First, suppose we allow for more than one infectious con- 
tact between individuals. In this case, compartmentalizing 
nodes by their number of unmatched stubs would be ill- 
advised. Indeed, when allowing for more than one infectious 
contact, the propagation dynamics cannot be made equiva- 
lent to a dynamical process performing on the fly neighbor 
assignment. Opting for a local description where nodes are 
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sorted according to the number and state of their neighbors, 
as in rf27l l28ll . would yield a more accurate description of 
the dynamics. Here, the basic state variables would read 
[^1 Y2]i s i,j s j,{t), where is and ij represent the number of neigh- 
bors on Pi that are respectively 1 -susceptible and 1 -infectious 
(same for j$ and ji regarding agent 2 on T 2 ). Note that track- 
ing the number of 1- and 2-recovered neighbors is optional. 

Second, the leaky partial immunity rule could be replaced 
by perfect partial immunity lfl8ll . In the latter, each g- 
susceptible node that is of g -state Y g regarding the other vi- 
ral agent g has a probability cr^ of being perfectly immune 
to agent g. This scenario could be analytically modeled by 
introducing a fourth state compartment for each viral agent, 
V g , denoting ^-susceptible nodes that are perfectly immune to 
agent g. 

Finally, our analytical approach could also be made com- 
patible with SIS dynamics after slight changes in the com- 
partmentalization scheme. Because SIS dynamics cannot be 
made equivalent to a dynamical process performing on the fly 
neighbor assignment, it would once again be necessary to opt 
for a local description where nodes are sorted according to the 
number and state of their neighbors 11271 12811 . 

IV. VALIDATION THROUGH CASE STUDIES 
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FIG. 3. (Color online) Total agent 1 incidence Ri(oo) versus the in- 
fectious contact rate B t for a single viral agent SIR epidemic in two 
host populations with different network topologies. Population A is 
characterized by a Poisson distributed network, while population B 
has a power-law distributed network. Both networks feature the same 
mean degree (ki) ^ 3.5. Note that near/Jj = 0.66 (vertical line), both 
epidemics have a similar incidence. Points and error bars correspond 
to the mean and standard deviation of Monte Carlo simulations; solid 
curves are the predictions computed from our analytical approach. 



use this particular/?i value throughout this section when com- 
paring together both host populations. Moreover, Fig.|3]shows 
that our mean-field approach is able to capture with great ac- 
curacy the behavior of the model in the single-agent case. 



In order to analyze the behavior of the model and to validate 
the accuracy of our mean-field approach, we consider two host 
populations, A and B, which are under the threat of an agent 
1 outbreak. Each host population is characterized by a dif- 
ferent network Fj . For population A, the network Fi is fairly 
homogeneous and features a Poisson degree distribution, 

Pi(h) = Ci^—-,0<k 1 <20, (18) 
k\ ! 

where A\ = 3.5 and is a normalization constant. We as- 
sume that the population B displays more heterogeneous con- 
tact patterns, and has a network Fj that follows a power-law 
degree distribution: 

Pf(/fci) = Cf/fcp , 1 < h < 20 , (19) 

where ti = 1.5 and is a normalization constant. The pa- 
rameters are chosen such that both networks have the same 
mean degree {k\) ^3.5 while featuring different level of het- 
erogeneity. 

Figure |3]shows the phase diagram for the total incidence of 
agent 1, defined by 

Y i,j,k 

as a function of the infectious contact rate B\ in the case where 
agent 1 propagates alone in the host populations A and B. 
R\{oo) is computed by setting e 2 = 0, in which case the dy- 
namics reduces to the case of a single viral agent SIR epi- 
demic. Note that for B\ = 0.66, both epidemics reach approx- 
imatively the same level of incidence. For this reason, we will 



A. Delayed intervention 

In response to an epidemic menace, intervention strategies 
involving the propagation of a second viral agent in the host 
population may be employed to control the outbreak of the 
undesirable agent. Examples of such strategies include the 
spread of preventive information in the context of an emerging 
disease, or the dissemination of countermeasures to minimize 
the damages of a computer virus outbreak. In this section, we 
use the framework of interacting epidemics on overlay net- 
works to analyze the efficiency of a hypothetical delayed agent 
2 intervention on the outbreak of agent 1 in the host popula- 
tions A and B. 

We consider a case of unidirectional immunity. We assume 
that agent 2 is not affected by agent 1, i.e., cr? = 1 for all 
X e {S,I,R}, while nodes that are either 2-infectious or 2- 
recovered benefit from a given level of immunity to agent 1, 
i.e., cr^ = 1 and o-[ — <x^ = cr. In both host populations, the 
network T 2 , on which agent 2 propagates, is characterized by 
a power-law degree distribution, 

P 2 (k 2 ) = C 2 k 2 T2 , 1 < k 2 < 40 , (21) 

where t 2 = 1 and C 2 is a normalization constant. For now, 
we assume that the overlap between the networks Fi and T 2 is 
random, and that there is no degree correlation between them, 
such that P(k) = P\(k x )P 2 {k 2 ). 

a. Full immunity. The total agent 1 incidence R\{°°) in 
host populations A and B as a function of the delay r between 
epidemics is illustrated on Fig. |4(a)| for the case of full (cr = 0) 
and partial (cr = 0.5) immunity. Let us consider the case of 
full immunity in the first place. For cr = 0, we observe in 
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FIG. 4. (Color online) (a) Total agent 1 incidence ^i(oo) in the host 
populations A and B versus the delay r before an agent 2 interven- 
tion (/? 2 = 1) providing full (cr = 0) or partial (cr = 0.5) immunity 
to agent 1. (b)-(c) Phase diagram showing Bi(oo) in both host pop- 
ulations as a function of the delay t and the infectious contact rate 
Bi of an agent 2 intervention providing full immunity to agent 1. 
B[ = 0.66 in all figures. In both population, agent 2 spreads on a 
power-law distributed Ti network that has a random overlap with the 
T[ network. Points and error bars correspond to the mean and stan- 
dard deviation of Monte Carlo simulations; curves are the predictions 
computed from our analytical approach. 



Fig. |4(a)| that R\(°°) increases much faster with t in popula- 
tion B than in population A. This means that one disposes 
of a much shorter time to react efficiently if the population 
in which agent 1 spreads features an heterogeneous structure, 
such as a power-law distributed network. Figures |4(b)| and 
|4(c)| which also include the dependency of Ri(oo) on the in- 
fectious contact rate B2, further confirm this observation. 

The results discussed in the last paragraph can be explained 
by looking at the time evolution of the epidemics. In Fig. [5] 
the prevalence at time t of agent 1 and 2, 



hit) = J] ^[/i^l^W . hit) = Y J Yj [Xlh]i * (t) 

Y i,j,k X i,j,k 



(22) 



is illustrated for increasing values of the delay t in the case of 
full immunity. For small values of t (t = and 1), agent 2 is 
able to inhibit the initial phase of the agent 1 epidemic in both 
host populations. At intermediate t values (r = 5), the agent 
1 epidemic is still strongly restrained by agent 2 in popula- 
tion A, while it has almost enough time to run its course com- 
pletely in population B. As r is further increased (r = 10), 
the effect of the intervention becomes minimal in both host 
populations. Figure [5] shows that the time scale of the agent 
1 epidemic in a given host population is crucial in determin- 
ing the efficiency of a delayed intervention. Since this time 
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FIG. 5. (Color online) Time evolution of the agent 1 and 2 preva- 
lence, /[(f) and hit), in host populations A and B for various values 
of the delay r before an agent 2 intervention providing full immunity 
(cr = 0) to agent 1. Infectious contact rates are B\ = 0.66 and 82 = 1. 
Points corresponds to the mean of Monte Carlo simulations; error 
bars are of order of magnitude of the symbol size. Solid curves are 
the predictions computed from our analytical approach. 



scale decreases with increasing network heterogeneity 112911 . 
this explains why much smaller values of t are required in 
population B to achieve an efficient intervention. 

b. Partial immunity. Let us now consider the case of par- 
tial immunity (cr = 0.5). For the host population A, Fig. |4(a)| 
shows that the behavior of R\i°°) versus t is quite similar to 
the case of full immunity. The total agent 1 incidence in- 
creases slightly faster with the delay, which can be attributed 
to the fact that the fraction of the population reached by agent 
2 before agent 1 is now partially vulnerable to the latter. How- 
ever, the picture is drastically different for the host population 
B. Even when r = 0, agent 2 is unable to inhibit the agent 1 
epidemic, which reaches half as many nodes as it would reach 
without any intervention. 

This phenomenon is also a consequence of network het- 
erogeneity. Consider a node of degree £1 on Fi that is still 
1 -susceptible and has acquired partial immunity to agent 1. 
If it is eventually involved in infectious contacts with all its 
neighbors on F\, the probability that it remains 1 -susceptible 
in the end is (1 - cr) kl , which decreases exponentially with 
k\. Therefore, leaky partial immunity has a weaker effect on 
high-degree nodes. This observation, together with the fact 
that the high-degree nodes are preferentially involved in the 
early phase of an epidemic If29ll . explains why agent 1 is able 
to invade a significant fraction of population B even when a 
short-delay intervention is attempted. 

Finally, Figs. |4(a)| and [5] show an excellent agreement be- 
tween the analytical predictions computed from our mean- 
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field approach and the outcome of Monte Carlo simulations of 
the dynamics. As we mentioned in Sec. IIII Al we expect our 
mean-field approach to be exact for configuration model over- 
lay networks of infinite size and random overlap. The small 
divergence between the predictions of our approach and the 
mean values computed over Monte Carlo simulations can be 
attributed to finite-size effects, such as stochastic extinction at 
early times and the restriction on self-loops and repeated links 
when generating the networks. 



B. Overlap and degree correlation 

In the previous section, we assumed that the overlap be- 
tween the networks H and ^ was random, and that there ex- 
isted no degree correlation between them. We now relax this 
assumption and look at the effect of overlap and degree corre- 
lation. 

We consider a scenario very similar to that of the previ- 
ous section, where agent 1 (disease) and agent 2 (interven- 
tion) propagate simultaneously in the host populations A and 
B. We assume that the intervention is instantaneous (r = 0) 
and grants full immunity to agent 1 (cr = 0). The degree dis- 
tributions of the network T\ for population A and B are respec- 
tively given by ( fT8l and < fT9b . Moreover, we now suppose in 
both populations that the degree distribution of Y2 is identical 
to that of H, i.e. P 2 (k) = P\(k) = p(k). 

a. Overlap versus degree correlation. In order to be able 
to isolate the respective effects of overlap and degree correla- 
tion, we distinguish between three different configurations: (i) 
random overlap and no degree correlation, in which case the 
system is built from P(k) = p{k\)p{k.2); (ii) random overlap 
and full degree correlation, in which case we build the sys- 
tem from P(k) = p{ki)6^,ki', (iii) full overlap, in which case 
p(K) = 5 Kl flS Kl flp(Kb) is used to generate the system. 

The bifurcation diagrams obtained for each host population 
using the three above configurations are illustrated in Fig. [6] 
Let us compare in the first place the two configurations where 
the overlap is random. In both host populations, we observe 
that agent 1 manages to invade the system more easily in the 
case where there is no degree correlation between the net- 
works. This is a consequence of the fact that the high-degree 
nodes of a network are more likely to be infected in an out- 
break. If the high-degree nodes on Fi and T2 correspond to the 
same individuals, both viral agents will preferentially com- 
pete for their infection. In the scenario considered here, it will 
therefore be harder for agent 1 to invade the system if the high- 
degree nodes on Fi are preferentially immunized by agent 2. 
The difference in heterogeneity between the host populations 
A and B explains why the effect of degree correlation is much 
stronger in B. 

Because full overlap implies full degree correlation, the 
fully overlapping configuration must be compared with the 
case of random overlap but full degree correlation in order to 
isolate the effect of overlap. In Fig. [6] we observe that a con- 
figuration with full overlap facilitates the epidemic of agent 1 
in the system. Since Fi = Y2 in the case of full overlap, nodes 
reached by agent 1 that are not already part of the agent 2 out- 




(a) Population A 




(b) Population B 

FIG. 6. (Color online) Total agent 1 incidence Ri(oo) versus the in- 
fectious contact rate B t in host populations A and B for different net- 
work configurations (■: random overlap, no degree correlation; 0: 
random overlap, full degree correlation; ♦: full overlap). In all cases, 
the agent 2 intervention (B2 = 1) is undelayed (t = 0) and provides 
full immunity (cr = 0) to agent 1 . Points and error bars correspond to 
the mean and standard deviation of Monte Carlo simulations; curves 
are the predictions computed from our analytical approach. 



break are more likely to be connected to other nodes that have 
not yet been reached by agent 2. For a given joint probabil- 
ity distribution P(k), this explains why invasion of agent 1 is 
easier in the fully overlapping case. 

Note that the magnitude of the difference between R 1(00) 
in all the different configurations of degree correlation and 
overlap is larger when the time scales of the agent 1 and 2 
epidemics are comparable. When B\ is too low, agent 1 can 
barely invade the system in all configurations, while when 
B\ » B2, agent 1 is almost unaffected by agent 2 because 
the former spreads significantly faster. 

Our results corroborate the previous findings of Funk and 
Jansen. In I20I1 . they considered the case where two processes, 
the first granting full immunity to the second, propagate sub- 
sequently on overlay networks. They showed that the epi- 
demic threshold of the second process increases with positive 
degree correlation between networks with random overlap, 
while it decreases with increasing overlap for a given joint de- 
gree distribution P{k). In our work, variations in the epidemic 
threshold are hard to analyze because both processes spread 
simultaneously in the host po pula tion. Our results generalizes 
the previous conclusions of ICOII to the case of dynamically 
interacting processes on overlay networks. 

b. Partial overlap. In Sec. IIII Bl some approximations 
were introduced in our mean-field approach to allow for arbi- 
trary overlap between Fi and To. To investigate their validity, 
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FIG. 7. Total agent 1 incidence R | (oo) in the host population A versus 
the percentage of overlap Q between Y { and T 2 . The parameters are 
a = 0, t = 0, and /3\ = fii = 1. Points and error bars correspond 
to the mean and standard deviation of 2500 Monte Carlo simulations 
performed with N = 100 000. The solid curve corresponds to the 
predictions of the analytical approach developed in Sec. HUB I The 
dashed curve corresponds to the analytical predictions obtained after 
making the substitutions given in Eqs. J24) . 



we now consider the general case where the system is built 
from the following joint degree distribution 

P(k) = p{K X + K b )d^h * - QT , (23) 

where Q is the mean percentage of overlap in the system. 

In Fig. |7l we compare the outcome of Monte Carlo simu- 
lations with the analytical predictions of the approach devel- 
oped in Sec. IIIIBI (solid curve) for increasing values of Q in 
the host population A. We see that our approach becomes less 
accurate as Q increases. This mainly comes from the fact that 
nodes can be granted additional unmatched stubs on j\ and j2 
after a neighbor assigment on yt>, which causes the dynamics 
on j\ and j2 to looses its exact character. Figure Q however 
shows that our approach becomes accurate again as Q gets 
very close to one. While this may appear counterintuitive at 
first, this is because an almost exact correspondence with the 
propagation dynamics is reobtained at Q = 1 . Since there are 
no links on the j\ and j2 networks at Q — 1, the i and j in- 
dices are initially zero for each node. Thus, they are only used 
to indicate the number of new stubs that are granted after in- 
fectious contacts on ji,, i.e., there is no mix between real and 
supposed unmatched stubs on j\ and j2- 

Finally, it is interesting to note that for the special case of 
full immunity, an approach that handles partial overlap more 
accurately can be obtained. Indeed, since agent 2 provides full 
immunity to agent 1, we may use the following substitutions: 



for Y e {S,I,R}. While this modified approach yields a better 
description of the dynamics for the epidemic scenario consid- 
ered here (dashed curve in Fig. |7), we emphasize the fact that 
it is not valid for the general case of an interaction rule involv- 
ing partial immunity. 

V. CONCLUSION 

In this paper, we have introduced a general model where 
two viral agents propagate simultaneously on overlay net- 
works and interact dynamically via a mechanism of leaky par- 
tial immunity. Exploiting a correspondence between prop- 
agation on networks and a dynamical process performing 
on the fly network generation ||2lll . we have developed a 
network-based compartmental formalism in which nodes are 
sorted at any time according to their state and number of un- 
matched stubs. Unlike previous work based on bond percola- 
tion lfl6l[T8l - t20ll . our analytical approach gives the complete 
time evolution of the system. By direct comparison with full 
Monte Carlo simulations of the model, we have demonstrated 
that it is able to capture with great accuracy the dynamical 
interaction between epidemics occurring simultaneously on 
overlay networks featuring various level of heterogeneity, de- 
gree correlation, and overlap. 

Our analytical approach is highly versatile and may be ap- 
plied to numerous scenarios of diverse nature. Here, we have 
considered a hypothetical delayed intervention scenario, in 
which an immunizing process is disseminated in a host popu- 
lation to hinder the propagation of an undesirable process. We 
have discussed the influence of the delay and the level of im- 
munity on the intervention efficiency in host populations fea- 
turing homogeneous and heterogeneous network structures. 
Moreover, we have shown that positive degree correlation in- 
creases the efficiency of the intervention, while overlap facili- 
tates the invasion of the undesirable process. 

Finally, this work highlights the power of the general mod- 
eling scheme presented in I2TI1 . from which our formalism 
stems. Part of ongoing research focuses on the application 
of these guidelines to include more realistic features in our 
model, such as community structured networks ll30ll . Other 
directions for future research include the application of our 
analytical approach to investigate specific case studies, such 
as the influence of other sexually transmitted infections on the 
spread of HIV & 
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S-I. NOTATION 



TABLE I. Notation used in this document 



P(kuh) 


Joint degree distribution of the system of overlay networks consisting of T t and T 2 . 


p(Ki,K 2 ,K b ) 


Joint degree distribution of the three non-overlapping networks decomposition consisting of y lt y 2 , and y b . 




Mean fraction of nodes in the host population that are of 1 -state X\, of 2-state Y%, have i unmatched stubs 
on 1*1, and j unmatched stubs on To at time t. 




Mean fraction of nodes in the host population that are of 1 -state X\, of 2-state Y 2 , have i unmatched stubs 
on y\ , j unmatched stubs on y 2 , and k unmatched stubs on y b at time t. 


g-, g 


The index g e {1,2) denotes the viral agent currently under consideration. The other viral agent is denoted 
by the index g € j2, 1}. 


atg 


Recovery rate for viral agent g. 




Infectious contact rate for viral agent g. 




Fraction of nodes in the host population that are initially g-infectious. 


< 


Probability of successful transmission of viral agent g after an infectious contact involving a g-susceptible 
node of g-state Fg regarding the other viral agent g. 




Probability of failed transmission of viral agent g after an infectious contact involving a g-susceptible node 
of g-state F| regarding the other viral agent g (a^ = 1 - crj). 


® g 


Probability that an unmatched stub chosen at random on network y g belongs to a g-infectious node. 


®b, g 


Probability that an unmatched stub chosen at random on network y b belongs to a g-infectious node. 


K s 


Probability that an unmatched stub chosen at random on network y b belongs to a g-infectious node whose 
g-state regarding the other viral agent g is such that a transmission of agent g may occur later between it 
and a node of g-state Yg, 


— Y 


Probability that an unmatched stub chosen at random on network y b belongs to a g-infectious node whose 
g-state regarding the other viral agent g is such that a transmission of agent g may never occur later between 
it and a node of g-state Yg (@ b g = & b , g - ® Y bg ). 


or 


Probability that an unmatched stub chosen at random on network y b belongs to a node whose g-state 
regarding the other viral agent g is such that a transmission of agent g may occur later between it and a 
node of g-state Y%. 


b,g 


Probability that an unmatched stub chosen at random on network y b belongs to a node whose g-state 
regarding the other viral agent g is such that a transmission of agent g may never occur later between it and 

a node of g-state Y g (<D^ = 1 - ). 



The following summation convention is used throughout this document: 

oo oo oo 

Z-Z- Z = Z- Z = Z- Z = Z- Z = Z 

X Xe{S,l,R] Y Ye{SJ,B] i i=0 j j=0 k k=0 
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S-II. OVERLAY NETWORKS WITH RANDOM OVERLAP 

Let 

©i = ; — and ©2 = : — . (26) 

For the case of overlay networks with random overlap, the following system of ordinary differential equations (ODEs) governs 
the time evolution of the [X\Y 2 ]ij compartments: 

4[SiS 2 ]y = A®i [0'+l)5f[5i5 2 ] (!+w - *[SiS 2 ]y] 

u? (27) 



-rlSihh = A©i [(» + l)^i[5i/ 2 ](i + i U - i[S i/ 2 ]y] 

w? (28) 

- a 2 [S 1 7 2 ] y +^ 2 ® 2 (J + l)af [5i5 2 ] i0+ i) +y 8 2 (l + © 2 ) [<J + 1)[S i/ 2 ]i U+ i) - y[S ihlj] 

J t \- S i R 2h =A®i f(i + l)af[5 1 /? 2 ] ((+1)i - i[Si*2]y] 

(29) 

+ a 2 [5 i/ 2 ]y + j8 2 © 2 [U + 1)[S 1^2]/ 0+ l) " j[S iRliij] 

ilhSih = - oriC/iSdy +A©i(* + [SiS 2 h+»j + + ®i) [(' + IH^fc+iy - *[/iS 2 ]y] 

^ (30) 

+ £ 2 02 [(; + l)^[/i5 2 ]/ 0+ l) - 7[/lS 2 ]y] 

^[/i/Jy = " ariCVdy + + l)^[Si/ 2 ]( i+ i)j +^i(l + ©i) [(i + l)[/i/ 2 ] ( , + iy - i[hh] U ] 

at (31) 

- a 2 [/i/ 2 ]y +y82©2(; + \)(j{[hS 1 ] im) +j8 2 (l + 2 ) [(j + l)[/i/ 2 ] i0+ i) - jlhhhj] 

-rUiRih = ~ axihRilj +A®i(* + l)<rf [Si*2](^i)7 + ®i) [(«' + l)[/i^2]( I+ i» - i[/i* 2 ]y] 

uf (32) 

+ a 2 [/i/ 2 ]y +^©2 [(j + l)[/i^2]; a+ i, - jUMj] 
^■[RiSih =ari[/iS2]y+A®i [(» + l)[*iS 2 ] (i+ i y - i[*iS 2 ]y] 

(33) 

+ /? 2 ®2 [U + [J?iS 2 ] i0+1) - j[J?iS 2 ]y] 
-[*i/ 2 ]ij = ari[/i/ 2 ]y +y8,©i [(i + l)[Rih\M)j - i[*i/ 2 ]y] 

\J'-rJ 



-[RiR 2 ]ij = adhRilj +/3i®i [(» + l)[^i*2]«+iw - i[*i*2]y] 

(35) 

+ a 2 [^i/ 2 ]y +yS 2 02 [(j + l)[RiR 2 ]tu+i) ~ ilWj] ■ 
This system of ODEs satisfy the conservation of nodes, namely 

The dynamics of agent 1 is initialized at time t — with the following initial conditions: 

[S iS 2 ]y(0) = (1 - ;) [/i5 2 ]y(0) = j) [RiS 2 ]ij(0) = 

[Si/ 2 ]y(0) = [/i/ 2 ]y(0) = [«i/ 2 ]y(0) = (37) 

[5 i/? 2 ]y(0) = [/ii? 2 ]y<0) = [fli/?2]y(0) = . 



The dynamics of agent 2 is initialized at time t — r > by making the following substitutions: 

[SiS 2 ]y(r) -» (1 - e)[SiS 2 ]y(T) ['iS 2 ]y(r) -> (1 - e2)[/i5 2 ]y(T) [^i5 2 ], 7 (t) -» (1 - ft)[/?i5 2 ]y(r) 
[S i/ 2 ]y(r) -> £2[5 i5 2 ]y(r) [/i/ 2 ]y(r) -» e 2 [/i5 2 ]y(r) [/?i/ 2 ]y(r) -> e 2 [/?i5 2 ]y(r) . 



(38) 



S-III. OVERLAY NETWORKS WITH ARBITRARY OVERLAP 



Let 



0, s 



®4,1 = 



YjYTji,j,ki[hY2]ijk 
YiX'J' Jji',f,k' i'lX'xY^li'j'k' 

'Ly'Li,i,kk[hY2\ijk 



09 = 



and ® b2 = 



lixliUjjttfXiIiiijk 



We further define the following probabilities: 



YiX',Y> k'lX^Y^li'fk' 
Hi,j,kk[IiS 2]ijk 







if 7 = 5 



if Y = I 



if Y = R 



W -0*! ifF = 5 



Zz',f Tji'jM k'lX^Y^li'fk' 

TjxZi,j,k k [XiS2]ijk 

T*X',Y> Tii'J'Jc' k'[X l Y 2 ]i'j'k' 



o 



if y = s 



if Y = I 



if y = /? 



— y 
and <I> fc , 



Zuyt(*ISi/2]</t+W2]ijit) 



Yji',j',k> k'lX^Y^li'j'k' 
Hi,jJik[S lh]ijk 



YjX'J' Tii'J'.k' k'\X x Y^\i'j'k' 





if y = 5 



if y = / 



if y = /? 



— F 

©4,2 = 







fc,2 



s 



© w -0i 2 







/>.2 



if Y = S 

if y = / 



©£ 2 ify = « 



'Y,xY,i,j,k{k[S l X 2 \ ijk +k[hX 2 \ ij k) 

TiX'.Y' Tii'J'Ji' k'lX^Y^li'j'k' 
~Lx Zi,j,k k \- S lXl]ijk 

TjX'J' ILi'jM k'lX^Y^li'fk' 



if y = 5 



if y = / 



if y = R 



— Y 

and <D i2 = 



i-o£ 2 ify = 5 

l-<t>£ 2 ify = 7 

i-of 2 ify = « 



14 



For overlay networks with arbitrary overlap, the following ODEs govern the time evolution of the [X1I2L7A compartments: 
^[SiS 2 ];./a- = A®i [(' + lW\[SiS2h+W ~ i[S\Si\ijk\ 



+#(*+ IF? (©i,i[5i5 2 ]/ -i)(* + i) +®l, l [SiS2]u(k + i))-Pi®b,ik[S l S 2 ]i 
+fa®2 [(J + l)^[SiS 2 ]iy + i)* -;'[Si5 2 ], 7 J 

+ /? 2 (fc + Dol (®l2[SiS 2 ] (i _ w+1) +0l 2 [5iS 2 ] ! ^ + i ) )- J 8 2 Aj2 fe[S 1 5 2 ] 



y* 



^[5i/ 2 ], 7 a =A®i [(' + i)^[5i/ 2 ] ( / + i,>t - HSihljk] 



- a 2 [S ihlijk + p2®2<J + [S iSilu+Dk +&(1 + ® 2 ) [(; + l)[S ih]i(M)k - j[Sih] tjk ] 
+ f3 2 (k+l)oj (®l 2 [S l S 2 ] (i -l)j(k+l) + ®b,2i S l S 2 ] ;,/(*+ 1) ) 



+ Pi(k+ 1) 



-\)j{k+\) + (^fc^ + ©fe,2)[ 5 lh]ij{k+\) 



■P2{l+®b,2)k[Slh}ijk 



-[5^2]^. =A®i [(* + [5i/? 2 ](i + i)/* - ^l^ly*] 



+ A(fc+l)^(®M[Sltf2]iO-l)(i+l)+®^^ 

+ a 2 [S ih] ijk +P2®2 [U + 1)[S iR 2 ]iu + m - W i R 2]ijk] 

+ p 2 {k + 1) (@l 2 [S iR 2 h-i)Kk + i) + ®l, 2 [S i^2]/ M+ i)) -P2®b,2k[S iR 2 ~]ijk 

-r[IlS2]ijk = ~ ai[/lS 2 ]^ +/?i©i(/ + l)(7i[SlS 2 ](M)jk + ®l) [U + WlS 2 ] {M)jk - iUlSrfijk] 



-fril + Qb^kUiSAijk 



{Ki + ®U)^ s ^u-i)(k + i) + (® b ,i + ® b ,i)UiS 2 ]m + i) 

+ /3 2 ® 2 [{J + l)o^[/i5 2 ] !(7+w - j'[/i5 2 ] i7 i] 

+ /3 2 (A: + IW ! 2 (©i^hS^y^ + ©fc i2 [/i5 2 ], 7((t+1) ) -p 2 ® b ,2k[hS 2 \ 
J^h\ijk = - axihhljk +Pi®i(i + l)cr[[S Ah+Djk +fi t (l + ®i) [(» + l)[V2](i + i);* - j[V2]yt] 
+ y3,(*+ (®ii[5i/ 2 ],o-i)(*+i) +®I,i[5i/2]y ( i + i)) 



ly* 



+ /3i(fe+l) 



- J 8l(l+®*,l)^ / l / 2] i7 i 



«1 + ®i,l)t/l/2](0-lXh-l) + +®6,l)[/l/2]y0h-l) 

- oraC/iJal/j* +^2® 2 (j + l^l/iSdim-i)* + ©2) [(i + - Ahhhjk] 

+ /3 2 (k+l)cr 2 (®' b2 [I l S 2 ] (i-l)yCt+l) + ©fa, 2 l7lS 2 ], 7 (£ + i) 



H,2 + ®k)[V 2 ]( ;-l) 7 (A+l) + (^,2 + ®fo,2)[ / l / 2]« 



;(*+!) 



-A(l+®M)W2]ij* 



^[/i^ 2 ], 7 a = - ari[/i*2]y* +A®iO' + iK^i^lf+Dji + + ®i) [(»' + l)[Ji*2]« + iy* - itWo*] 



+ 



O'-lXA+l) + ®i,l[^ 1^2]i 7 (A+l) 
— R —R 



) 



Kl +KlY I l R ^U-W+D + (%,! +® b , l )lhR2]ij(k + l) 

adhlriijk +/3 2 ® 2 [(j + l)[hR 2 ] iU+l)k - j[hR 2 ]ij k ] 

j3 2 (&+l)(®£ 2 [/l*2](i -l)M+D + ®fe, 2 [A^2], 7 (A+l) -^2®A,2*[A^2]y* 



-/3i(l +06,i)*[/i/?2]yt 



(44) 



(45) 



(46) 



(47) 



(48) 



(49) 



(50) 
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-fi[RiSdijk = aAhSiljk + /3i©i [(» + l)[RiS 2 ] ( M)j k - *[*iS 2 ]y*] 

+ p 2 ® 2 [u + iW2[RiS 2 ] iU+ m - ARiSzhjk] 

+ /3 2 (k + l)o* (0^ 2 [^i5 2 ](,-i U(it+ i) + ©t.^i^yxit+i)) - A®w*I*iS 2 ]i>* 
J t l R ih]ijk = adhhljk + /3i©i [(i + l)[Rihh + i)jk - mhljk] 



+ /3 2 (k+l) 

17?, if,l, a = 



K,2 + ®m)[*i^](M)XH-D + K.2 + ©y^i^yci+D -AO + ©*, 2 W 2 ]y* 



-j-[*i*2]y* = aitWy* [(» + UMdenvt - mRiljk] 



+fii&+ l)(©?,i[^2]io-ix-t+i) + ®M[*i*2]y(* + i)) -A®W*[*i*2]ii* 

+ Or 2 [«l/2]ij« + /3 2 ©2 [U + l)[^1^2]i 0+ l)-t - /[*1^2]^] 

-l)Kk+l) + ®6,2[^1^2]y(i+l) I -/?2®ft,2^[^1^2]/;* 



Once again, it is possible to verify that this system of ODEs satisfy the conservation of nodes, 



(51) 



l)(®l,i[Rih]iu-w+i) + ©w[*i'2]y(* + i)) -fiiQb.iti.Rihlijk 

- a 2 [RiI 2 ] ijk + /3 2 ® 2 (j + VaflRiSzliU+W +j8 2 (l + ® 2 ) [(j + l)[RihMj+W ~ RRihljk] (52) 
+ /3 2 (k+ (0^[*iS 2 ]«-il»+i) +©i ;2 [^i5 2 ] y(*+i) 



(53) 



x,r u,k 

The dynamics of agent 1 is initialized at time t — with the following initial conditions: 

[5 iS 2 ]y*(0) = (1 - eOfKi, j, k) lhS 2 ] ijk (0) = e lP (i, j, k) [*iS 2 ]y*(0) = 

[Sih] ijk (0) = UihhkiO) = [«i/ 2 ]yi(0) = (55) 

[S iR 2 ] U k(0) = [/i*2]y*(0) = [*i*2]yt(°) = . 

The dynamics of agent 2 is initialized at time t — t > by making the following substitutions: 

IS iS 2 ] ijk ( T ) -» (1 - e 2 )[S iS 2 ]</t(T) [/iS 2 ]^(t) -> (1 - (&)[/iS 2 ]y*(T) [*iS 2 ]</t(T) -> (1 - £ 2 )[RiS 2 ] ijk (T) 

IS Jiljkir) -» e 2 [5 iS 2 ] ijk (T) [hhljkir) -» ^DOTyiM [*i/ 2 ]y*(T) -» e^i^y^r) . po; 



